ANISOTROPIC FINITE ELEMENTS WITH HIGH ASPECT RATIO 
FOR AN ASYMPTOTIC PRESERVING METHOD FOR HIGHLY 
ANISOTROPIC ELLIPTIC EQUATIONS 
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Abstract. The concern of this work is the generalization of an Asymptotic Preserving method for 
the highly anisotropic elliptic equations presented in |14| . The limitations of the method introduced 
there in are omitted by the introduction of a stabilization term in the Asymptotic Reformulation. 
Furthermore, anisotropic error indicators and mesh adaptation algorithms are proposed and tested 
allowing to reduce considerably the number of mesh points required to achieve prescribed precision. 
Reported meshes have maximum aspect ratio greater than 500. 
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1. Introduction. Anisotropic problems are common in mathematical modeling 
of physical problems. They appear in various fields of application, such as flows in 
porous media '17] , semiconductor modeling pT| , quasi-neutral plasma simulations 
[TT] . image processing [28, 29_, atmospheric or oceanic flows [57] and so on, the list 
being not exhaustive. The direct motivation of this work is related to numerical sim- 
ulations of strongly magnetized plasma such as internal fusion plasma of tokamak 
13 , atmospheric plasma [121 HD] or plasma thrusters |1 . In this context a strong 
magnetic field is defining the anisotropy direction. Fast rotation of charged particles 
around magnetic field lines is causing a large number of collisions in the plane per- 
pendicular to the magnetic field. On the other hand the motion in the direction of 
the field is rather undisturbed. In consequence the particle mobility depends on the 
direction and may differ by several orders of magnitude. Anisotropy ratio 1/e can be 
as high as lO^''. 

The main difficulty associated with these anisotropic problems is that they are 
singular in the limit e — > 0. On the discrete level this is manifested by very bad 
conditioning of linear systems obtained by a direct discretization of the problem for 
e <C 1. In this paper we propose an approach based on the Asymptotic Preserving 
reformulation introduced initially by Shi Jin in |18j . Our approach is an extension of 
the method proposed in a previous paper |14) to the case of more general anisotropy 
field structure (such as closed field lines). 

The model problem we are interested in, reads 

-V • A,Vu^ = f infl, 

n-A^V-u^^O onrjv, (1.1) 
u'^ — on r^i , 

where C is a bounded domain with boundary dil — r^urjv and outward normal 
n. The direction of the anisotropy is given by a vector field i?, where we suppose 
divi? = and B 0. The direction of B shall be denoted by the unit vector field 
b = B/\B\. The domain boundary is decomposed into Tjj := {x G dfl \ b{x) -71 = 0} 



Universite de Toulouse, UPS, Institut de Mathematiques de Toulouse, F-31062 Toulouse, France 

1 



and Tjsf :— dfl\TD. The anisotropic diffusion matrix is then given by 

A, = -A^^b®b+{Id-b(E)b)Aj_{Id-b(g)b) . (1.2) 

The scalar field j4|| > and the symmetric positive definite matrix field Aj_ are of 
order one while the parameter < £ < 1 can be very small, provoking thus the high 
anisotropy of the problem. The system becomes ill posed if we consider the formal 
limit e — > 0. It is thus very ill conditioned for e ^ 1. 

This problem has been studied before in the Asymptotic Preserving context. A 
special case of anisotropy direction aligned with one of the coordinate axis was ad- 
dressed in |12j . A generalization of this approach was presented in , where the 
problem with curvilinear anisotropy field was reduced to one with the anisotropy di- 
rection aligned with the coordinate system by a change of variables. Another work 
pU] proposed a different generalization based rather on the introduction of Lagrange 
multipliers. This resulted in a considerably bigger linear system but allowed to avoid 
a necessity of change of variables which could be troublesome for time dependent 
anisotropy direction. Finally, a different method presented in 14] allowed to reduce 
considerably computational cost without any adaptation of the coordinate system. All 
those methods shared the same drawback: they didn't allow more complex geometries 
such as the presence of closed field lines. 

In this paper we introduce yet another Asymptotic Preserving scheme, improving 
the idea presented in |14j and removing the restrictions on the anisotropy direction by 
a simple penalty stabilization technique. Furthermore, the anisotropic error indicator 
is presented and the mesh adaptation algorithm developed in order to optimize the 
number of mesh points required to obtain a prescribed error. 

The outline of the paper is following. Section [2] contains a definition of the 
problem and introduces the Asymptotic Preserving reformulation. Section|3]describes 
an anisotropic error indicator and mesh adaptation algorithm. They are both tested 
and the numerical results are provided. 

2. Problem definition. We consider a two dimensional anisotropic problem, 
given on a regular, bounded domain il C M^, with boundary dil. The direction of the 
anisotropy is defined by the vector field b(x), which satisfies the following hypothesis 

Hypothesis A The field b{x) is derived from a vector field B{x) ~ \B{x)\b{x), 
satisfying div B{x) = and \b{x)\ = 1 for all a; G $7. Moreover, we suppose that 

b € [c°^[n)Y. 

Given this vector field 6, one can decompose now vectors w G M^, gradients V(/), with 
4>{x) a scalar function, and divergences V • w, with v{x) a vector field, into a part 
parallel to the anisotropy direction and a part perpendicular to it. These parts are 
defined as follows : 

v\\ := {v ■ b)b , w_L := {Id — b® b)v , such that v = v\\ + I'd. , 

V||0 := (6 • V(/))6, Vj_(t>:= iId-b<^b)V(j}, such that V (/) = V + V , 

V|| • u := V • W|| , Vj_ ■ V := V ■ vj_ , such that V • t; = V|| ■ v + Vj_ ■ v , 

(2.1) 

where we denoted by (g) the vector tensor product. With these notations we can now 
introduce the mathematical problem, the so-called Singular Perturbation problem, 
whose numerical resolution is the main concern of this paper. 
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2.1. The Singular Perturbation problem (P-model). The objective of this 
paper is to introduce an efficient scheme for the precise (e-independent) resolution of 
the following Singular Perturbation problem 

-iV|| • (A||V||0^) - • (A^Vi^^) - / inn, 

(P) { i?i||-(A||V||0^)+nx-(A^Vi0^) = O ondn,,,Udn,.,t, (2.2) 

= on drio , 

where n is the outward normal to il and the boundaries are defined by 

dnD = {xedn\b{x)-n = 0}, (2.3) 
dn,^ ^{xedn\ b{x) • n < 0}, (2.4) 
dnout ^{xC,dVL\ b{x) • n > 0}. (2.5) 

The parameter < e < 1 can be very small and is responsible for the high anisotropy 
of the problem. We shall assume in the rest of this paper the following hypothesis on 
the diffusion and source terms 

o 

Hypothesis B Let f G L'^{Q,) and diljj ^ 0. Furthermore, the dijfusion coefficients 
A|| e L°°{n.) and A± e M.dy.d{L°°{n)) are supposed to satisfy 

< Ao < A||(a;) < ^1 , fa.axen, (2.6) 

AjLix)b{x) ^ A]_{x)b{x) ^ , fa.axeQ, (2.7) 

MM'^ < v*'Aj_{x)v < , yveR'^ with v ■ b{x) = and f.a.a x e n. 

(2.8) 

As we conceive to use the finite element method for the numerical resolution of the 
P-problem, let us put (2.2 1 under variational form. For this let V be the Hilbert space 

V := {(^ e H\n) I <i,\Q^^ = 0} , (0, V)v := (V,i0, V,i^)l2 + e(y ^<i>, ViV)^^ . 

We are searching thus for 0^ e V, solution of 

a,|(0^^)+eai(</)^7^)=£(/,^), VV-eV, (2.9) 

where (., .) stands for the standard L'^ scalar product and the continuous, bilinear 
forms a|| : V x V M and aj^ : V x V M are given by 

a||(</),7/;) / A,, V||0 • V,, V rfx , aj^(0» := / (A^lV^l-/-) • V^lV' f^a^ . (2.10) 

Thanks to Hypothesis B and the Lax-Milgram theorem, the problem (2.2) admits a 
unique solution (f)^ £V for all fixed e > 0. However, the numerical resolution of (2.2 1 
is very inadequate for e <C 1. When e tends to zero, the problem reduces to 

-V|| • (yl||V||(j!)) = in 17, 

nil . (A|| V||<j!)) = on 8^,^ U d^ouu (2.11) 

0° = O OTidVtD- 
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This is an ill-posed problem as it has an infinite number of solutions 4> ^ S, where 



g = {(j,eV\W\\cf> = 0}, (2.12) 

is the Hilbert space of functions, which are constant along the field lines of b. On the 
discrete level this is manifested by a very bad conditioning of the system for small 
values of e. However, as shown in [ID], the solution (j)'^ £ V converges to 0° G t/, a 
unique solution of 



(L) [ Aj_S/±(j)° ■\7±^|Jdx ^ [ fipdx , \f^p e G . (2.13) 
Jn Jn 

2.2. The Asymptotic Preserving approach (AP-model). Let us introduce 
a so called AP-formulation, which is a reformulation of the Singular Perturbation 



problem (2.2), permitting a "continuous" transition from the (P)-problem (2.2) to the 



(L)-problem (2.13), as e 0. The AP-formulation was introduced and is a subject of 
more detailed analysis in a separate publication [T3] . We will shortly recall the results 
of the previous studies. For this, each function shall be decomposed into two parts: 
constant part along the anisotropy direction and a part containing fluctuations. The 
constant part converges to the limit solution and the fluctuating to as e (see 
also tMj). 

Let us introduce the following Hilbert space: 

A:={qe L^{n) /V\\q G L^in) and q\an,^ = 0} (2.14) 
(g,u))^ = (Viiq, Viiw) , yq,weA. (2.15) 



Let 0^ be a solution to the Singular Perturbation problem (2.2) and set 



+ eq'^ with p"^ £ Q and q"^ G A. This decomposition is unique and we observe 

( a_L(p'',w) + ea_L(gr^,w) + a||(g^,w) = (/,w) Vw G V , 



(2.16) 

aii{p'^ ,ui) — VwG^, 



or equivalently 



r a{r,v) + (1 - e)a,| (g^ v) = (/, v) G V , 
(AP) (2.17) 
[ a|| (0'^, w) = ea|| (q^, w) Ww^A, 

with the bilinear form a{v, w) deflned as 

a{v,w) = / AVw • Vw. (2.18) 
Jn 

The matrix A is given by 

A = ® 6 + (/d - 6 (g) h)A^{Id -h®h), (2.19) 

and is e independent, A = Ai. 

The above formulation is the Asymptotic Preserving reformulation based on the 
Micro Macro decomposition. 
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2.3. The stabilized Asymptotic Preserving approach (AP-model). The 

Asymptotic Preserving approach presented above has some hmitations originating in 
the choice of the vector space A. Note that in the previous paper the uniqueness of 
was ensured by setting g*^ to on the Ti„ boundary under hypothesis that every 
field hue of b has its beginning on Tin and an end on Tout- In other words, more 
complex geometries, like for example closed field lines are not permitted. In this 
paper we propose a new way of providing the uniqueness of q^ which overcomes the 
limitations of our previous method. The idea is based on the penalty stabilization 
method introduced in [5] for the Stokes problem. 

Let us propose a new Asymptotic Preserving method: find {(j)^ ,q^) G V x V such 

that 

r a(0^ i;) + (1 - £)a|| (g^ v) = (/, v) Vt; G V , 

{APS) { (2.20) 
{ a|| (,^^ w) = eaii (<7^ w) + Y.Ker^ hj, AVq' -Vw Vu; G V , 

where hfc denotes the size of the element K. Note that now, instead of seeking q^ E A 
we are looking for G V. Existence and uniqueness of the above problem can be 
easily proved by the Lax-Milgram theorem. 

3. Numerical method. This section concerns the discretization of the Asymp- 



{APS) 



h 



totic Preserving formulation ( 2.20 1, based on a finite element method. The anisotropic 
error indicator is introduced and the obtained numerical results are studied. 

Let us denote by V/t C V and Ah C A the finite dimensional approximation 
spaces, constructed by means of Pi finite elements. We are thus looking for a discrete 
solution G Vh X Ah of the following system 

a{(t>l,Vh) + (1 - e)a\\{ql,Vh) = {fh,Vh) ^Vh G Vh , 

a\\i<l>l:Wh) = £a\\iqh,Wh) + J2KeT^^K iK^'^lh '"^^h ^Wh G Vh ■ 

(3.1) 

3.1. Adaptive finite elements with large aspect ratio. We now propose an 
adaptive finite element algorithm. The goal is to build successive triangulations with 
large aspect ratio such that the relative estimated error of the function (j)'^ — p^+eq'^ in 
the H^{n) norm is close to a preset tolerance TOL. For this purpose, we introduce an 
error indicator which requires some further notations. This error indicator measures 
the error of the numerical solution (j)^ in the directions of maximum and minimum 
stretching of the triangle. The goal of the adaptive algorithm is then to equidistribute 
the error indicator in the directions of maximum and minimum stretching, and to align 
the directions of maximum and minimum stretching with the directions of maximum 
and minimum error. We refer to [24l [23j El [HI [16] for theoretical justifications. 

For any triangle K of the mesh, let : K —i' K he the affine transformation 
which maps the reference triangle K into K. Let Mk be the Jacobian of Tk that is 



X = Tk(x) = MKX + t 



K- 



Since Mk is invertible, it admits a singular value decomposition Mk — R]^AkPk, 
where Rk and Pk are orthogonal and where Ak is diagonal with positive entries. In 
the following we set 



Ak = ( ^^Q^ 1 and Rk 



l.K 
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with the choice Xi k > ^2,k- A simple example of such a transformation is xi = Hxi, 
X2 — hx2, with H > h, thus 



M, 



K 



H 
h 



Ai 



K 



H, A 



2.K 



h, ri,K 



r2,K 



see Figure |3.1[ In other words Ti k and r2.K are the directions of maximum and 
minimum stretching, while Xi_k and X2,k measure the amplitude of stretching. 

Let Ifi ■ HQ{fl) Vh he a Clement or Scott-Zhang like interpolation operator. 
We now recall some interpolation results due to [151 HH 123 • 

Proposition 3.1. There is a constant C = C{K) such that for all v e H^{fl), 
for all K € Th, for all edges e of K , we have 



\v - Ihv\\L2(n) < C [Xlj^iri^KG k{v)vi^k) + >^I^k{'^2,kG k{v)v2,k)) 



1/2 



(3.2) 



\v - IhV\ 



L2(e 



< Ch 



1/2 / Xi^K 



K 



X 



{ri,KGK{v)ri^K) + 'l^{r2,KGK{v)v2,K) 



^2,K 



' X2 



Ai, 



K 



1/2 



1/2 



V(w - Ihv)\\L^(K) < C I j^{ri^KGKiv)ri^K) + (r2,xGi^ (w)r2,7<) 



(3.3) 



(3.4) 



Here h^ — diam K, Xi^K o^nd Vi^k o.tc given by (3.1), and Gk{v) denotes the 2x2 
matrix defined as 



Gk{v) 



( 



K 



dv 
dxi 



dx 



K 



dv 
dxi 



dv 

0X2 



dx 



\ 



K 



dv 
dxi 



dv 
dx2 



dx 



K 



dv 
dx2 



dx 



(3.5) 



Proof. The first estimate is in Proposition 3.1 of [T^, the second estimate is in 
Proposition 2.2 of [16], the third estimate is in Proposition 2.5 of [22]. □ 

The results of Proposition [3T] are now used to derive an anisotropic error indicator 
for the Asymptotic Preserving reformulation. The error is first related to the equation 
residual. The Clement interpolant is introduced. Then the anisotropic interpolation 
results are used. Finally, a Zienkiewicz-Zhu error estimator is used to approach the 
error gradient. 

Let e = 0'^ — 0| and — ^ The following error estimate for the Asymptotic 
Preserving reformulation (2.171 holds. 

Proposition 3.2. There exist a constant C depending only on the interpolation 
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constants from Proposition \3.1\ and not on the mesh size nor aspect ratio such that 
/ AVe • Ve + (1 - e)£ / A|| Vye, • V|[eq + (1 - e) h\ j AVe, • Ve, 

Jn Jn Kern 

C E (ll/ + V-(AV<^D + (l-£)V|| •(A||V||gD||, 



\L^{K) 



1 l-e 



1/2 ML" ' Y-ft ■-iwi.'yaf.) • ^^1/2 

1/2 



ZA2 _^ z/\2._ft: 



X (A?j^(ri,A'Gif(e)ri,K) + \l,K{^2,KGK{e)v2,K)) 
+ (1 - e) ( l|V|| • (^||V||(0f, - £9D)llL^(i^) + -^||[A||V||g^ • n]\\L^aK) 

3/2 I 



X {>AA^i.KGK{eq)vi,K) + A2 ;^(r2,KGA(e,)r2,A))'/' . (3.6) 

Here [■] denotes the jump of the bracketed quantity across an internal edge, [■] = for 
an edge on the boundary dflo, [■] is set to twice the imposed flux on the dQinU dflout 
and n is the unit edge normal in arbitrary direction. 

Proof. Setting v — e in the AP reformulation ( 2.17[ ) yields 

a(e,e) + (1 - e)a||(e,eq) = (/, e) - a(0^,e) - (1 - e)a\\{ql,e). (3.7) 

Now, since ay (0^ — eq'^ , Cq) = X^iferh ^i" !k ^^^^ ' we obtain 

a||(e,eg) = £a||(e,,e,) + V / AVg'' • Ve, - a,[ (0^ - £g^, e,) (3.8) 
and hence 

/ AVe • Ve + (1 - e)£ / A|| Vye, • Vye^ + h\ \ AVe, • Ve, = 

/ fe- f AV0^Ve -(l-e)/ A,, V||<7, • Vpe +(!-£)/ A||V||(0^, - egfj • Vpe, 

-{l-e)Y,hl [ kVql ■ Ve,. (3.9) 

For any w € we have 

(/, v) - a{(pl,v) - {1 ~ e)a\\{ql,v) 

= (/,"- Ihv) - a(0^, V - Ihv) - (1 - e)a|| (qf,, -y - hv) 

= E ( / (/ + V-(AVr,) + (l-£)V|| •(A,iV,igD)(t'-^^^') 

+ 1 f [KV<t>l-n]{v-hv) + ^ [ [A\\V^^qi-n]{v-Ii,v)\ . (3.10) 



Furthermore, for any w G A the following holds true : 



= <^\\i<Ph -£qh,w- Ihw) - h\ I ■V{w- Ihw) 



-hj^ V • (AVg|)(w - h^v) + hj, / {AVql ■ n){w - hw) . (3.11) 

Jk JdK J 

Now, choosing v — e, w = Cq and using the Cauchy-Schwartz inequality together with 
the interpolation results of the Proposition 3.1 the following is obtained: 

/ AVe • Ve + (1 e)£ / A|| V||e, • V||eg 
Jn Jn 

+{l-e) hi f AVe,.Ve,<C^ ( ||/+V.(AV^^J + (l-£)V||.(^l|V||gDllL^W 
+ 1 f x y \\[AVcj^ln]\\^.^oK) + ^ ( . \ )'^'||[A||V||g|.n]|U.o^)) 



1/2 



1/2 



/ 1 / h \ 

+ (l-£) i ||V,i-(yl,iV,i((/)|-£g^))||L2(A) + 2 ( aikIa J 

X (A?,K(i-i,AGK(e,)ri,K) + A|^(r2,A-GK(eg)r2,A))'^' (3.12) 
where C = C{k). Since J^^ ^|| V||e^ • V||eq > and 

^i.Khj^ <hK < \2,Khf^, (3.13) 

the inequality (3.6) holds true. □ 

Remark 4. Note that the above result does not contain any terms inversely 
proportional to e as it involves matrix A rather than A^. The standard anisotropic 
error indicator for an anisotropic diffusion problem studied in ]23l \25j takes form: 

f A,Ve-Ve<G^ j ||/ + y • (A,V</)|) + -^||[A,V0^ n]|U.(a^) ) 

X {\lK(ji,KGK{e)vi^K) + A^,^(r2,AGK(e)r2,A))'^' , (4.1) 

thus it involves terms of the order - . While this error indicator remains valid it is of 
no practical use for small values of e. Indeed, the remeshing algorithm which aims in 



keeping the error indicator close to a given value would yield meshes with mesh size 
proportional to e. 

Remark 5. In the case of e ~ I the above error indicator reduces to the standard 
anisotropic error indicator for a diffusion problem studied in [] : 



f AVe • Ve < C V ( 



||/ + V-(AV0^J 



2A, 



1/2 
2.K 



X i>^l,K{^i,KGK{e}ri^K) + A^_K(r2,KGK(e)r2,if))^^^ . (5.1) 



Estimate (3.6) is not a usual a posteriori error estimate as it involves 

<Ph and q 



and q"^ 



the right hand side. If we can guess <j)^ — 0| and q^ — qf^, (3.6 1 can be used to derive an 
anisotropic error indicator. In order to do that, we introduce an error estimator based 
on the superconvergent gradient recovery, namely Zienkiewicz Zhu like error estimator 
O 1301 131] in its simplest form as defined in [51 [55], i.e. the difference between V(j)l^ 
resp. V^l and an approximate projection of V(/)| resp. V(7| onto : 



( 



{I - n,, 



dxi J 



(5.2) 



where II/j is the projection operator which builds values at vertices P from constant 
values on triangles using the formula 



(P) 



E 1^1 



tria. K 



E 

tria. K 

peK 

E 

tria. K 

, PGK 



K\ 



\K\ 



dxi 
dx2 



\K 



\K 



Z-Z like error estimator is asymptotically exact for a parallel meshes and smooth 
solutions 12 US]. Our error indicator is obtained by replacing the matrices Gi<-(e) and 
Gif(eg) by approximate ones GK{4>h) ^nd GKiqfi) defined by 



G 



{vf^i^Dfdx 

K 

ri?^mril^{4>V)dx 



vrmvrmdx 



K 



K 



ivrmrdx 



(5.3) 



K 



The anisotropic error indicator defined on each triangle K takes the form 
[ri1c{'PUh))f = (^11/ + V • (AV<^|) + (1 - £)V|| • (A||V||(zD|U2(K) 

2-^2, if / 
X (A?,^(ri,KGK(<^|)ri,if) + A2^(r2,KGK('ADr2,x))'^^ 

+ (1 - e) ( l|V|| • (^||V||(0^ - £gD)IU^(K) + — I7^II[^I|V||(<^| - eql) ■ n] 1^2(3^) 
+ A2,if l|V • (AVg^)||i2(;^) + \l[l\\kWql ■ nWmaK^ 

1 /2 

X (Atif(ri,KGif(g|)ri,K) + A^_^(r2,ifGK(?Dr2,if)) • (5.4) 

Introducing 

P<I>,K = 11/ + V • (AV<^D + (1 - £)V|| • V||g^)|U2(^) (5.5) 
2A2 2A2^ 

(5.6) 

1/2 



{viK{<Pl Qh))) = P<P,K (xlK{Yl,KGK{4'l)ri,K) + Xl,K{Y2,KGK{<Ph)^2,K) 



(5.7) 



and 



Pa,K = (!-£) (^||V|| • (^||V||(</.| - eqlML^K) (5.8) 
+ 1 1 [^11 ^11 (-^h - ■ n]\\L-(dK) (5.9) 

+ A|;f ||V • (AV5|)|U2(^) + X\%\\AVql ■ n\\ L-^(dK) , (5.10) 

/ n2 / _ ^ sl/2 

yitMh Qh))) = Pa,K i^XlK{ri,KGK{(l>h)n,K) + XI,k{->^2,kG K{4'h>2,K)) 

(5.11) 

allows to introduce a more compact notation 



{r,i{<l>Ul))f = {viM<Pllh))y + {<K{<l>Uh))f ■ (5.12) 

5.1. Adaptive algorithm. The goal of our adaptive algorithm is to build a 
triangulation such that the error is cquidistributed in the direction of the maximal 
and minimal stretching of triangles and the relative global error indicator is closed to 
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prescribed tolerance TOL. We have 

^ , ^ - < 1 



0.75 TOL < \ ^ < 1.25 TOL. (5.13) 



with 

2 



(77^(<^|,g|)) = ^ (ry^|(0|,gD) ■ (5.14) 



tria. K 



A sufficient condition to satisfy ( 5.13 1 is to build a triangulation with large aspect 
ratio such that 



OJb^TOL'^ ,9 1.252TOr 



7VT 



for each triangle K, where NT is the number of triangles of the mesh. Since the mesh 
generator BL2D mesh generator used in our simulations f7^ requires data on the mesh 
vertices rather than on the triangles, we need to translate the above local triangle 
condition into a condition for mesh for the mesh vertices. Let us introduce a point 
defined error indicator : 

4 



ri^i^l,€)= E W^l^ll)) (5.15) 



tria. K 

peK 



and hence 



E('?p(<^^''?^))' = 3E(^k(<^L9D)'. (5.16) 



P K 

Therefore, the following local condition holds 



^0.75'TOL'J^ iVr.P < {v^{c^f.,cil)y < ^1.25'tOL' (5.17) 

where NV is a number of mesh vertices. Then, we define vfpi'f'h' ifi)' with i = 1,2 
at the mesh nodes 

{vtpi'l^llDy = E (^''A' {pl,KGKm+PlKGK{ql)) V^.k) ■ (5.18) 



tria. K 



The value of vfpi'Ph'^h) represents the error in the direction of the maximum and 
minimum stretching of the triangle K. We note that the point error indicator is 
bounded by 



{viA^lqf.))' + {vtpi'l^lqf.))') ■ (5.19) 



< 2 
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The mesh adaptation algorithm can be summarized as follows. For all vertices 
P of the mesh Vi'pi't'hT'lh) '?2^p('^/i' 9|) '^^'-^ computed. Furthermore, we compute 
Ai,p and A2,p as an average of the Ai^^ and X2,k of the neighboring triangles K. 

The input data for the BL2D mesh generator is computed: the stretching amplitude 
hi.p, i = 1,2 and the direction of the anisotropy 6p. In the first step new hi^p are 
obtained. For every mesh point P, if 

^{vtpm,€)y < j^0.75'TOL^ \Vcj>l\'^ (5.20) 
then hi^p is set to 3/2Ai^p. If 

^{vtpi^hQDy > ^j^l.25'TOL^ (/^ imi')' (5-21) 

then hi^p is set to 2/3Ai^p. Otherwise, hi^p is set to Xi^p- 

In the second step of the mesh adaptation the new anisotropy direction is found. 
For every mesh point average matrices Gp{<f)\) and Gp{qf^) are calculated. The angle 
is set to the angle between the eigenvector corresponding to the largest eigenvalue 
of the matrix 

pIkGk{<I>1)+pIkGk{cQ (5.22) 

and the Ox direction. Finally, new mesh is generated using the BL2D mesh generator. 

5.2. Simplified error indicator. The anisotropic error indicator introduced in 
the previous sections involves the term (g^). This means that the perpendicular 
derivatives of g^ will play role in the error estimation procedure. This is not necessarily 
desirable since in some cases this may result in mesh over-refinement in the direction 
perpendicular to the anisotropy direction. That is to say the adaptive algorithm could 
continue to refine the mesh in the perpendicular direction without any increase of 
precision in This is why we propose an alternative approach where the simplified 
error indicator is related only to the residue of the first equation and the matrix 

GKm- 

{rilc\<t>l,ql)y = (^\f + V • (AVc/.^) + (1 - £)V|| • (A,, V||g^)|U.(K) 

X {xlK(v^,KGKWh)^i,K) + X\Av2,KGKmV2,K)\ ''^ , (5.23) 



or in more compact notation: 

(^i^('^L'ZD))' = {viK{^lql))y- (5.24) 
As in the previous section the nodal simplified error indicator is defined: 

E('?^^('^'^'9^))' = ^J2{^K^i^l'€)y, (5.25) 

P K 

(»?^^p(<^LaD)'= E PlK>^lKri,KGKmri,K. (5.26) 



tria. K 
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The obtained adaptive algorithm is almost the same as before. Only now 'qfp is 
replaced by a simpHfied version rj^p, the coarsening criterion is slightly changed : if 

'^{<p{^Ul)f < j^0.75^TOL^ iV^^J^y (5.27) 
then hi^p is set to 3/2Ai,p. If 

then p is set to 2/3Aj,p. Otherwise, hi,p is set to A;. p. Finally, the mesh anisotropy 
direction is aligned with the largest eigenvalue of the matrix Gk{4>%)- 

5.3. Numerical results. 

5.3.1. Numerical study of the effectivity index and the convergence of 
the stabilized AP scheme. Let us define 

1/2 



\ 1/2 



the Z-Z error estimator, the anisotropic error estimator and the simplified error indi- 
cator. We also define 

^zz 

ei^^ = —L , (5.32) 

' (la -^Ve ■ Ve + t(l - s) J„ A,V,,, ^ V,e,)"/2 ' P'^^' 

the effectivity indices. 

We test the robustness of the error indicators and the convergence of the stabilized 
AP scheme in the following test case. Let f2 = (0, 1) x (0, 1), the anisotropy direction 
is given by 

6 = — B=( "^^^ ~ C0S{WX) +TT \ ,^ 

\B\' \ 7ra(y2 - y) sin(7rx) J ' ^ ' ' 

Note that we have B ^ in the computational domain. The parameter a describes 
the variations of the anisotropy direction. For a = the anisotropy is aligned in the 
direction of x coordinate. Wo set A± = ^|| = 1. Now, we choose 4>^ to be a function 
that converges to the limit solution (j)^ as e ^ 0: 

(jp = sin (ttj/ + a(j/^ — y) cos(7ra:;)) , (5.36) 

(j)^ = sin (ttj/ + a{y^ — y) cos(7rx)) + e cos (2itx) sin (ttj/) . (5.37) 
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Finally, the force term is calculated accordingly, i.e. 



/ = -Vi • (AiVi0^) - -V|| • (A||V||0^). 

We study the effectivity indices on the unstructured meshes for constant and 
variable anisotropy direction (a = and a = 2 respectively) and for small and large 
anisotropy (e = 1 and e = 10"""^° respectively). 



hi-h2 








l|V(</'|-r)llL^(0)/l|V(/)||U.(^,) 


0.1-0.1 


1.05 


2.53 


2.53 


1.5 X 10-1 


0.05-0.05 


1.02 


2.54 


2.54 


7.7 X 10-2 


0.025 - 0.025 


1.01 


2.54 


2.54 


3.9 X 10-2 


0.0125-0.0125 


1.00 


2.53 


2.53 


1.9 X 10-2 


0.00625 - 0.00625 


1.00 


2.53 


2.53 


9.8 X 10-^ 


a = 0, e = 1 


hi-h2 




ei^ 




||V(</.|-r)l|L^(0)/||V(/)||U.(f,) 


0.1-0.1 


0.99 


4.74 


4.68 


8.0 X 10-2 


0.05-0.05 


0.99 


4.78 


4.71 


4.1 X 10-2 


0.025 - 0.025 


0.96 


4.76 


4.67 


2.1 X 10-2 


0.0125-0.0125 


0.93 


4.89 


4.65 


1.1 X 10-2 


0.00625 - 0.00625 


0.87 


5.08 


4.68 


6.1 X 10-^ 


a = 0, e 10-1° 


hi-h2 








l|V(<^|-r)llL^(0)/l|VrJlL=(n) 


0.1-0.1 


1.05 


2.54 


2.54 


1.5 X 10-1 


0.05 - 0.05 


1.02 


2.54 


2.54 


7.7 X 10-2 


0.025 - 0.025 


1.01 


2.54 


2.54 


3.9 X 10-2 


0.0125-0.0125 


1.00 


2.53 


2.53 


1.9 X 10-2 


0.00625 - 0.00625 


1.00 


2.53 


2.53 


9.9 X 10-3 


a = 2, e = 1 


hi-h2 




ei^ 




l|v(<^|-</>^)IU^(o)/||V0||U.(a) 


0.1-0.1 


0.99 


4.07 


3.99 


1.1 X 10-1 


0.05-0.05 


0.98 


4.24 


4.07 


5.4 X 10-2 


0.025 - 0.025 


0.97 


4.30 


4.09 


2.7 X 10-2 


0.0125-0.0125 


0.94 


4.41 


4.08 


1.4 X 10-2 


0.00625 - 0.00625 


0.90 


4.69 


4.19 


7.5 X 10-3 



a = 2, e = 10-1° 
Table 5.1 
Effectivity indices for isotropic meshes 



Table 5.1 shows the numerical results for isotropic unstructured meshes in different 
regimes. In the case of no anisotropy (e = 1) the Zienkiewicz-Zhu error estimator 
converges to true error as h goes to zero. The simplified and full effectivity indexes 
are the same and converge also to a constant value. In the case of small anisotropy 
(e = 10-1°) the effectivity index for Zienkiewicz-Zhu error estimator is close to one 
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for all testes isotropic meshes in the case of variable anisotropy direction. However, 
its value seems to decrease with the mesh size meaning that the estimator slightly 
underestimate the true error for fine meshes. The divergence is observed for a constant 
direction of anisotropy and small value of s. This shows that the Zienkiewicz-Zhu 
error indicator is not always equivalent to the true error. The stabilized Asymptotic 
Preserving scheme converges to the exact solution in all four cases with the optimal 
convergence rate. 

Table |5.2| presents the numerical results corresponding to the of large anisotropy 
aligned with the coordinate system. This time we are interested in the behavior of 
the error indicators when the mesh refinement is anisotropic. In the first table the 
mesh is refined in the direction perpendicular to the anisotropy direction with aspect 
ration ranging from 10 to 1280. In this case the Zienkiewicz-Zhu remains constant 
and close to 1. The relative error converges until the aspect ratio of 80 is reached. The 
effectivity index for the full error indicator increases from 6.28 to 15.6 with the mesh 
size until the aspect ratio reaches the value of 160. At the same time the effectivity 
index for the simplified error indicator is between 5.57 and 6.64. This suggests that 
the latter could perform better in the anisotropic mesh refinement. Its effectivity 
index does not seem to depend on the aspect ration wen the mesh is refined in the 
"right" direction (perpendicular to the anisotropy). 

Next, the influence of the mesh refinement in the "wrong" (parallel to the aniso- 
tropy) direction is performed. For aspect ration ranging from 1 to 16 the divergence 
of the ei^^ and the relative error is clearly observed. In fact, all effectivity indexes 
approach zero with the refinement. The last table displays the results of the con- 
vergence of ei^^ in the case of anisotropic mesh with aspect ration 4 and triangles 
aligned in the "wrong" direction. In this case, when the mesh is refined in both di- 
rection, the effectivity index for Zienkiewicz-Zhu error estimator approaches 1. The 
effectivity indexes of both error indicator diverge. 

5.3.2. Mesh adaptation. We now apply our adaptive algorithm to build a 
sequence triangulations in the following way starting from an isotropic unstructured 
grid with h — 0.02. At every iteration of the algorithm the error indicator is used 
to construct a subsequent mesh. We compare results of the simplified and full error 
indicators in various regimes: small and large anisotropy, b direction constant and 
variable. We focus on the resulting mesh size and error in the 7?^-norm as well as on 
the error convergence in terms of prescribed tolerance TOL. 



Let n = (0, 1) X (0, 1), the anisotropy direction is given by (5.351. We set A± = 
A|| — 1. We choose cj)^ to be a function that converges to the limit solution (fP as 



(f) — ii'm [Try -\- a{y ~ y) cos{7tx)) ^ ^ , (5.38) 



, Q , / '7Ty-\-a(y —y) cos^TTx) — 0.5\2 

(f)'^ ^ sin [ny + a{y — y) cos(7ra;)) e^*- ' ' -I- e cos (27ra;) sin (tt?/) . 

(5.39) 

Finally, the force term is calculated accordingly. The limit solution is nothing else 
than the limit solution from previous section multiplied by a Gaussian following the 
anisotropy direction. The parameter S controls the width of the exponential part. 
Setting S — 0.1 in our simulations yields a solution which has a strong gradient 
in the direction perpendicular to the anisotropy direction in a small subregion of a 
computational domain. The adaptive algorithm should be able to capture this strong 
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hi ii2 


ei 


ei 


ei 


II V((/>| - r)||L2(o)/||V?i|||i2(f^) 


0.1 - 0.01 


0.98 


6.28 


5.83 


7.8 X 10"^ 


0.1 - 0.005 


0.97 


7.68 


6.09 


4.3 X 10-3 


0.1 - 0.0025 


0.95 


10.6 


6.20 


2.5 X 10-3 


0.1 - 0.00125 


0.93 


14.2 


6.64 


1.7 X 10-3 


0.1 - 0.000625 


0.95 


15.6 


5.57 


1.6 X 10-3 


0.1 - 0.0003125 


0.98 


9.88 


3.94 


2.3 X 10-3 


0.1 - 0.00015625 


0.97 


9.20 


4.01 


2.2 X 10-3 


0.1 - 0.000078125 


0.98 


5.09 


2.95 


3.1 X 10-3 



Aspect ratio from 1:10 to 1:1280 



hi h2 








||V(</>|-(/.-)|U.(o)/||V(/<^IU.(n) 


0.1 - 0.1 


0.99 


4.74 


4.68 


8.0 X 10-^ 


0.05-0.1 


0.91 


4.56 


4.33 


1.2 X 10-1 


0.025 - 0.1 


0.32 


4.75 


4.04 


4.4 X 10-1 


0.0125-0.1 


0.005 


0.44 


0.36 


5.2 X IQi 


0.00625-0.1 


0.0002 


0.075 


0.059 


1.9 X 103 



Aspect ratio from 1:1 to 16:1 



hi-h2 








||V(</.|-r)llL^(o)/||V</.||U.(^,) 


0.025 - 0.1 


0.32 


4.75 


4.04 


4.4 X 10-1 


0.0125-0.05 


0.40 


6.66 


5.66 


2.0 X 10-1 


0.00625 - 0.025 


0.53 


8.86 


7.53 


7.4 X 10-^ 


0.003125 - 0.0125 


0.55 


9.54 


8.09 


3.6 X 10-2 


0.0015625- 0.00625 


0.69 


12.87 


10.09 


1.4 X 10-2 



Aspect ratio 4:1 
Table 5.2 

Effectivity indices for anisotropic meshes for a = and e = 10~^' 



variation of the solution and produce a mesh that is much finer in this subregion than 
in the remaining part of the domain. 

Small anisotropy e = 1, constant and variable direction of 6 (a = and 

a = 2) 

In the first two test cases the adaptive algorithm is studied in the e = 1 regime, 
i.e. when no anisotropy is present. In this case the two error indicators : full and 
simphfied are equivalent. 

Tables [STSl and [5^ show the results for b field with constant and variable direction 
respectively. The values in the tables are given after 15 iterations of mesh adaptation 
algorithm. In both cases the optimal convergence is obtained. The true error is clearly 
related to the prescribed error tolerance TOL and the node number is multiplied by 
4 every time TOL is divided by 2. The Z-Z effectivity index converges to 1 with TOL 
and the values of indexes for error indicators remain almost constant. This is not 
surprising since in this case the proposed error indicators reduce to the standard a 
posteriori error indicator studied before. The adapted meshes are presented on Figure 
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TOL 


err 


NV 




ei^ 




0.25 


0.096 


698 


1.03 


2.55 


2.55 


0.125 


0.048 


2457 


1.01 


2.54 


2.54 


0.0625 


0.024 


8834 


1.00 


2.57 


2.57 


0.03125 


0.012 


34587 


1.00 


2.54 


2.54 



Table 5.3 

error (err ), number of nodes (NV ) and effectivity indices for mesh iteration 15 for constant 
direction of b and e = 1 



TOL 


err 


NV 








0.25 


0.094 


785 


1.03 


2.58 


2.58 


0.125 


0.047 


2696 


1.01 


2.59 


2.59 


0.0625 


0.024 


10141 


1.00 


2.59 


2.59 


0.03125 


0.012 


39035 


1.00 


2.58 


2.58 



Table 5.4 

Relative error (err), number of nodes (NV) and effectivity indices for mesh iteration 15 
for variable direction of b and e = 1 



Numerical relative error obtained on the isotropic uniform mesh with h — 0.00625 
(31325 mesh points) give the relative error equal to 0.021, which is comparable with 
the results obtained for TOL — 0.0625. The adapted giving the same numerical 
precision are three times smaller. 

constant direction of & (a = 0), large anisotropy e = 10"^" 

In the next test case we consider large anisotropy e — 10~^° and aligned b di- 
rection. The simplified error indicator and the full error indicator are no longer 
equivalent. The results presented in Table [sTs] display the true error and effectivity 
indexes obtained by applying those two different algorithms. In this particular case we 
display results after 30 mesh adaptations. The number is bigger than in previous case 
in order to allow the algorithm to fully converge and exploit the reduced dimension- 
ality of this particular test. Note that in both cases the true error is comparable and 
converges with TOL. The Zienkiewicz-Zhu effectivity index is close to 1 for both error 
indicator. The aspect ratio for the smallest TOL studied is over 500. The simplified 
error indicator seems to perform better : the mesh size for the smallest TOL tested 
is three times smaller than for the full error indicator. The relative error is also 
slightly smaller for the simplified error indicator. The adapted meshes are presented 
on Figure [5^2] 

Numerical relative error obtained on the isotropic uniform mesh with h — 0.00625 
(31325 mesh points) give the relative error equal to 0.035, which is comparable with 
the results obtained for TOL = 0.0625. The adapted giving the same numerical 
precision are 115 (40) times smaller for the simplified (full) error indicator. 

variable direction of 6 (a = 2), large anisotropy e = 10^^° 
In the last studied test case we have applied the mesh adaptation algorithm to 
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(a) solution for constant b 



(b) adapted mesh 





Figure 5.1. Exact solution and meshes obtained after 15 iterations for TOL = 0.125 with e = 1 
for constant and variable direction ofb. 



the problem with large anisotropy with variable direction. Table 5.6 shows obtained 
results of numerical simulations. The simplified error indicator performs more effi- 
ciently than the full error indicator. Poor performance of the full error indicator for 
the smallest tolerance is caused by the perpendicular derivatives of which cause 
the over refinement in the direction perpendicular to the anisotropy direction. The 
resulting mesh is almost eight times bigger. For smaller values of the tolerance the 
difference in mesh sizes is much smaller and the meshes constructed for the full error 
indicator give slightly better precision. In both cases the Z-Z error estimator is close 
to 1. The adapted meshes are presented on Figure [5731 

Numerical relative error obtained on the isotropic uniform mesh with h = 0.00625 
(31325 mesh points) give the relative error equal to 0.038, which is comparable with 
the results obtained for TOL = 0.0625. The adapted giving the same numerical 
precision are 20 (10) times smaller for the simplified (full) error indicator. 
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(a) solution 



(b) full error indicator 



(c) simplified error indicator 



Figure 5.2. Exact solution and meshes obtained after 30 iterations for TOL = 0.125 with 
e = lO-""^" and constant direction ofb. 



TOL 


err 


NV 


(!l2\ 

\ Jmax 


(!l2) 


ei^^ 


ei^ 


0.25 


0.072 


272 


67 


12 


1.01 


3.20 


0.125 


0.037 


758 


86 


14 


1.01 


3.26 


0.0625 


0.018 


2435 


91 


17 


0.99 


3.32 


0.03125 


0.0093 


6642 


296 


23 


0.98 


3.28 



full error indicator 



TOL 


err 


NV 


(h:2\ 






ei^^ 


0.25 


0.060 


105 


130 


35 


1.01 


3.65 


0.125 


0.031 


271 


224 


48 


1.00 


3.49 


0.0625 


0.016 


652 


501 


87 


0.98 


3.74 


0.03125 


0.0076 


2018 


536 


106 


0.99 


4.07 



simplified error indicator 
Table 5.5 

Relative error (err), number of nodes (NV), maximum aspect ration ((^)max), average 
aspect ration ((j^)avg) and effectivity indices for mesh iteration 30 for constant direction ofb and 
e = 10-10 



6. Conclusion. A stabilized Asymptotic Preserving method for strongly aniso- 
tropic Laplace equation has been proposed and tested numerically. The error indica- 
tors including first order derivatives has been developed for this reformulated problem. 
Numerical experiments show the performance of the remeshing routine. The resulting 
meshes are considerably smaller by the factor from 3 to 115 than the isotropic uniform 
grids giving the same precision. The biggest gain is obtained for strong anisotropy in 
the constant direction. 
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